program celest3d
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
      use objects_com_M
      use controllo_M
      use celeste3d
      use cplot_com_M
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------

!-----------------------------------------------
!
!
!------------------------------------------------------------------
!               main program
!------------------------------------------------------------------
!
!     Initializations. Some variables are set to zeros and others are set according to the
!     input file.
!
      call initialization
!
!     main loop
!
      if (ncyc2 > 0) then
!
!
!
!l    >>>>>>>>>>>>>>> ncyc-loop <<<<<<<<<<<<<<<
!
         do ncyc = ncyc1, ncyc2
         !
            nh = mod(ncyc,nhst)
            nh = max(1,min(nhst,nh))
            t = t + dt
            write (*, *) 'TIME', ncyc, t, ncyc2
            told = t2
            !call second (t2)
            xx = (t2 - told)*rbjbkb
            if (ncyc == ncyc1) then
               costcc(nh) = 0.
            else
               costcc(nh) = (t2 - t2old)/(ibar*jbar*kbar)
            endif
            t2old = t2
!
!     ****************************************************
!
!     prepare for the next computation cycle
!
!     ****************************************************
!

            if (.not.testpar .and. (ncyc>1 .or. restart)) then
!
               call vinit
!
            endif
!

!			call Faraday
!

			if (eandm) then
			    !call maxwell_risolver
			    call faraday
			end if
!
            call budget
!
            if (lpr /= 0) then
!
               if (ncyc - 1==0 .or. mod(ncyc,50)-1==0) pflag = 1.
               if (pflag > 0.) write (kpt, 9000)
               write (kpt, 9010) ncyc, t, tmass, tmomx, tmomy, tmomz, racc, &
                  zacc, tie, tke, tbe, tle, numitv, numit, dvbmax, xx
               pflag = -1.
!
               call output
!
            endif
!     <<<<<<<<<<<<<<<<<<<<<< Solve Maxwell's Equations >>>>>>>>>>>>>>>>>>>>>>>>>
!
!
!     zero working arrays
!

            wate(:itdim,:5) = 0.0
            phipot(:itdim) = 0.0
!
if (fields) then
               error = 1.E-3
               relax = 1.
               tiny = 0.E-2

!
!	enforce div(E(t=n))=rho(t=n)
!
!            call maxwell_risolver
!            !
!            phibc = 0.
!            wk = 1.


!               call vtoctov (nvtx, ijkvtx, iwid, jwid, kwid, ibp2, jbp2, kbp2, &
!                  0, sdlt, cdlt, wkl, wkr, wkf, wke, phibc, ex, ey, ez, a11, &
!                  a12, a13, periodicx)
!            call faraday
    	if(divergence_cleaning) then
	       call div_cleaning
	    end if
!
!	Applying BC to the new E0
!
               call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, &
                  sdlt, cdlt, ex0, ey0, ez0, periodicx)
!
!	Smmothing the new E0
!
               phibc = 0.
               call vtoctov (nvtx, ijkvtx, iwid, jwid, kwid, ibp2, jbp2, kbp2, &
                  nsmooth, sdlt, cdlt, wkl, wkr, wkf, wke, phibc, ex0, ey0, ez0, &
                  a11, a12, a13, periodicx)
!
!
!    Electrostatic  potential (for graphics output, not needed for the cycle)
!
!
!
     if(iout(4)/= 0) then
               do n = 1, ncells
                  ijk = ijkcell(n)
                  wate(ijk,10) = fourpi*qdnc0(ijk)
               end do
!
               dumdt = 0.0d0
               dumscal = -1.d0
               wk = 0.
               phibc = 0.d0
               call Poisson(ncells, ijkcell, nvtxkm, nvtx, ijkvtx, itdim, iwid&
                  , jwid, kwid, precon, tiny, relax, ibp1, jbp1, kbp1, cdlt, &
                  sdlt, strait, dx, dy, dz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y&
                  , c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z&
                  , c6z, c7z, c8z, vol, vvol, wate(1,28), wate(1,1), wate(1,2)&
                  , wate(1,7), wate(1,8), itmax, error, wate(1,10), dumdt, &
                  dumscal, wk, phibc, ex, ey, ez, wate(1,4), wate(1,5), wate(1,&
                  6), phipot, wate(1,9), wkl, wkr, wkf, wke, periodicx)
!
               potavp(:ibp2*jbp2*kbp2) = potavp(:ibp2*jbp2*kbp2)*.99d0&
                    + phipot(:ibp2*jbp2*kbp2)*0.01d0
     end if
!
!
            if (eandm) then

            call maxwell_solver

            endif
!
            phibc = 0.
            wk = 1.
!
!
               call vtoctov (nvtx, ijkvtx, iwid, jwid, kwid, ibp2, jbp2, kbp2, &
                  0, sdlt, cdlt, wkl, wkr, wkf, wke, phibc, ex, ey, ez, a11, &
                  a12, a13, periodicx)
!
!
!
!     **************************************************
!
!	add an external electric field or other force
!	mapped into an effective electric field
!
!     **************************************************
!
!	call external_field

endif
!
!
!
!     **************************************************
!
!     solve the particle equations of motion
!
!     **************************************************
!
!
            if (testpar) then
!
               ncellsp = 1
               i = int(pxi(nphist+1))
               j = int(peta(nphist+1))
               k = int(pzta(nphist+1))
!
               ijkctmp(1) = (k - 1)*kwid + (j - 1)*jwid + (i - 1)*iwid + 1
!
            else
!
               ncellsp = ncells
!
            endif
!
            if (anydrift) then
!
               where = 1.
!      call counter(ncells,ijkcell,iphead,link,number,npart,where)
               call parmovgc
               where = 1.1
!      call counter(ncells,ijkcell,iphead,link,number,npart,where)
            endif
!
            if (explicit) then
!
!      call parmove
               if(relativity) then
               call parmovi_rel
               else
               call parmovi
               end if
!
            else

               if(relativity) then
               call parmovi_rel
               else
               call parmovi
               end if

!
!
!     particle control
!
     call control(sok,sko,cok,cko)
     write(*,*)'splitting: fail=',sko,'    succ=',sok
     write(*,*)'coalescen: fail=',cko,'    succ=',cok


            endif
!
!
            call ParInject
!
            if (testpar) call celindex (2, ibp1, 2, jbp1, 2, kbp1, iwid, jwid, &
               kwid, ijkcell, ijkctmp, ncells, ijkvtx, nvtxkm, nvtx)
!
!     **************************************************
!
!     interpolate particle quantities to the grid
!
!     **************************************************
!
            if (.not.testpar) then
               where = 5.
!      call counter(ncells,ijkcell,iphead,link,number,npart,where)
               call parcelc
!
!
               transpos = 1.
!
               call parcelv (ncells, ijkcell, nvtx, ijkvtx, nsp, iwid, jwid, &
                  kwid, ibar, jbar, kbar, dt, ijkctmp, anydrift, eandm, itdim, &
                  iphead, iphd2, link, transpos, qom, ico, qpar, pxi, peta, &
                  pzta, up, vp, wp, pmu, bxv, byv, bzv, wate, uptilde, vptilde&
                  , wptilde, mask, pixx, pixy, pixz, piyy, piyz, pizz, qdnv, &
                  number, j12x, j12y, j12z, jmag, jxtilde, jytilde, jztilde, jxs, &
                  jys, jzs,  pixysp2, piyzsp2, pixzsp2, pixxsp2, piyysp2, pizzsp2, &
                  pixysp1, piyzsp1, pixzsp1, pixxsp1, piyysp1, pizzsp1, &
				  uxsp1,uysp1,uzsp1,densp1,uxsp2,uysp2,uzsp2,densp2,avg)
!
!
               if (nrg_obj > 0) then
!
                  call parcelv_obj (ncells, ijkcell, nvtx, ijkvtx, nrg_obj, &
                     iwid, jwid, kwid, ijkctmp, itdim, iphead_obj, iphd2, &
                     link_obj, ico_obj, massp_obj, pxi_obj, peta_obj, pzta_obj&
                     , chip_obj, wate, chi, mv)
!
               endif
!
!
!
            endif
!
!

            !call second (t3)
            t3=0d0
!            This line below used to stop the simulation when a certain total simulation
!            time was reached (actual cpu wallclock)
!            if (t>=twfin .or. t3-stim>=tlm) go to 1010

!
            if (t2 - t1 < dcttd) then
!
               if (t + 1.E-10 < pttd) cycle
!
               pttd = pttd + dpttd
!
            endif
!l       (write all necessary data on tape (dump))
            t1 = t2
!
         end do
!
      endif
!
!l    >>>>>>>>>>>>>>> end of ncyc-loop <<<<<<<<<<<<<<<
!
      ncyc = ncyc2
!
!ll   4. terminate run:
!        -------------
      call output
!
 1010 continue
  212 format(i5,9(1p,e11.3))
!
!
      call finito
!
!
      call plotclose
!
      call restartout
!
      write (*, *) 'Im stopping from celest3d after ncyc=', ncyc

!
 9000 format('  cyc',8x,'t',6x,'mass',4x,'x-mom',3x,'y-mom',3x,'z-mom',6x,&
         'racc',6x,'zacc',5x,'dei',5x,'dek',5x,'deb',5x,'etotl',' itv',' it',2x&
         ,'dvbmax',4x,2x,'grinds')
 9010 format(1x,i4,1p,e9.2,1e11.3,1x,3e8.1,2e10.3,3e8.1,e10.3,i4,i3,e8.1,4x,e&
         8.1)
!
	stop
!
      end program celest3d

